Shock waves in strongly interacting Fermi gas 
from time-dependent density functional calculations 
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Motivated by a recent experiment [Phys. Rev. Lett. 106, 150401 (2011)] we simulate the collision 
between two clouds of cold Fermi gas at unitarity conditions by using an extended Thomas-Fermi 
density functional. At variance with the current interpretation of the experiments, where the role 
of viscosity is emphasized, we find that a quantitative agreement with the experimental observation 
of the dynamics of the cloud collisions is obtained within our superfluid effective hydrodynamics 
approach, where density variations during the collision are controlled by a purely dispersive quantum 
gradient term. We also suggest different initial conditions where dispersive density ripples can be 
detected with the available experimental spatial resolution. 
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A Fermi gas of atoms at unitary conditions, i.e. when 
the s-wave scattering length diverges [l|, is predicted to 
obey universal hydrodynamics, where the shear viscosity 
and other transport coefficients are universal functions of 
density and temperature. Experiments on the expansion 
of rotating, strongly interacting Fermi gas in the normal 
fluid regime reveal extremely low viscosity hydrodynam- 
ics 0, thus making a strongly interacting Fermi gas in 
the normal fluid regime a nearly perfect fluid with al- 
most frictionless mass flow. Remarkably, the superfluid 
behaves in rotating clouds experiments almost identically 
to the normal fluid one. Consistently with such observa- 
tions, theory predicts [3] that the unitary Fermi gas has 
zero bulk viscosity in the normal phase, whereas in the 
superfluid phase two of the three bulk viscosities vanish 

as. 

Nonlinear evolution of trapped cold gases is often char- 
acterized by the appearence of wave-like localized distor- 
sion of the trapped gas cloud. In addition to the well- 
known sound waves, shock waves often appear, which are 
characterized by an abrupt change in the density of the 
medium Shock waves are ubiquitous and have been 

studied in many different physical systems [(| 0] • Very re- 
cently the observation of nonlinear hydrodynamic waves 
has been reported in the collision between two strongly 
interacting Fermi gas clouds of 6 Li atoms at unitarity 
Q. During the collision, which has been experimentally 
imaged in real time @, when regions of high density 
move with faster local velocity than region of low den- 
sity, large density gradients develop when the two clouds 
merge with each other. Unlike the case of collision be- 
tween two initially separated BEC in a cigar shaped trap, 
where the shock waves were interpreted as purely disper- 
sive 0, the authors of Ref. Q invoke dissipative forces 
to avoid the " gradient catastrophe" in the unitary Fermi 
system. To describe their experimental data, they use 
an effective one-dimensional (ID) model based on time- 
dependent non-linear hydrodynamical equations, with a 
phenomenological kinematic viscosity term added to de- 
scribe dissipative effects, whose strength is used as a fit- 



ting parameter. A quite high value of the viscosity co- 
efficient is necessary in the description of the collision 
in order to reproduce the experimental data. The ori- 
gin of such viscosity, however, is not clarified in Ref. 
In fact the role of dissipation itself in cold Fermi gas 
at low temperature is questionable, since the ultracold 
unitary Fermi gas is known as an example of an almost 
perfect fluid, as discussed above. In this paper we offer 
an alternative explanation, where the collision process 
is dominated by dispersion effects Evidence sup- 

porting the notion that shock waves in ultracold Fermi 
gas at unitarity are dispersive rather than dissipative can 
be found in the simulations reported in Ref. Ill], based 
on the self-consistent solution of coupled Bogoliubov-de 
Gennes equations derived from the the zero-temperature 
time-dependent superfluid local density approximation 
12]. Our calculations are based on a single-orbital den- 
sity functional (DF) approach to the properties of uni- 
tary Fermi gas at low temperatures, which has been 
used recently to successfully address a number of prop- 
erties of such system }13l - ll8l ] . The advantage in using a 
single-orbital DF approach is that systems with a very 
large number of particles can be treated using only a 
single function of the coordinate, i.e. the particle den- 
sity. Thus, it might represent a viable alternative to the 
much more computationally expensive approaches (like, 
e.g., the Bogoliubov-de Gennes method) often used to 
describe superfluid Fermions. This is especially true in 
the case of 3D geometries, like the ones investigated here, 
where the BdG method would be prohibitively costly. 

In our extended Thomas-Fermi (ETF) density func- 
tional approach [lU [l4| the total energy of the unitary 
Fermi gas is given by 



E[n] = J d 3 r £(n, Vn) 
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Here n(r) is the fermions number density and U(r) is the 
confining external potential. The total energy functional 
£ contains a term proportional to the kinetic energy of a 
uniform non interacting fermions gas^plus a gradient cor- 
rection of the von-Weizsacker form [19| . In recent papers 
(13H15I . H3] we have determined the parameters £ and A 
by fitting Monte Carlo results [2(], HH for the energy of 
fermions confined in a spherical harmonic trap close to 
unitary conditions. The main conclusion of that work is 
that the values 



£ = 0.40 and A = 1/4 



(3) 



fit quite well Monte Carlo data of the unitary Fermi gas. 
In particular the chosen value for £ almost coincides with 
the experimental determination of Ref. [22j . 

The density functional ([IJ describes various static and 
dynamical properties of the unitary Fermi gas trapped by 
an external potential. The gradient term in the previous 
equation is found to be crucial to describe accurately the 
surface effects of the system, in particular in systems with 
a small number of atoms, where the Thomas-Fermi (local 
density) approximation fails [l3|, EE] ■ As a key result of 
the present paper, we show that when fast dynamical 
processes occur and/or when shock waves come into play 
such term is necessary also in the large N limit. 

The extended non-viscous and irrotational hydrody- 
namics equations deriving from the functional ([1} arc 
given by 



| + V.(„v) = 0, 



dv 2 d£ d£ 
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(4) 
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where n(r, t) is the time-dependent scalar density field 
and v(r, t) the time-dependent vector velocity field. If 
A = 0, then Eqs. (j4j and © reproduce the equations of 
superfluid hydrodynamics [23| by construction. 

Notice that equations (g]) and can equivalently be 
written in terms of a superfluid time-dependent nonlin- 
ear Schrddinger equation (NLSE) involving a complex or- 
der parameter [13| . We will numerically solve this NLSE 
equation to obtain the long-time dynamics of the col- 
lision between two initially separated Fermi clouds, as 
described in the following. Our goal is to simulate the 
experiments of Ref. [1] . We have used the Runge-Kutta- 
Gill fourth-order method 24j to propagate in time the 
solutions of the NLSE. To accurately compute the spatial 
derivatives appearing in the NLSE, we used a 13-point 
finite-difference formula [25j . 

Since the confining potential used in the experiments 
is cigar-shaped, we have exploited the resulting cylindri- 
cal symmetry of the system by representing the solution 
of our NLSE on a 2-dimensional (r, z) grid (of 500x2500 
uniformly spaced points). Of course, this choice, which 
greatly reduces the computational cost of the simula- 
tions, would not be able to describe possible, azimuthal 
dependent, transverse instabilities and vortex formation 



(like those observed in the collision between BEC clouds 
[3] ) . Although such features are apparently not observed 
in the collision experiments of Ref. (8j , we cannot rule out 
the possibility of transverse instabilities and vortex for- 
mation in Fermi gas clouds for different initial conditions 
than those investigated here. To describe the details of 
such structures, a full 3-D simulation is needed. 

In our simulations we tried to reproduce as closely as 
possible the experimental conditions of Ref. Q, which 
we summarize briefly in the following. A 50:50 mix- 
ture of the two lowest hyperfine states of 6 Li is confined 
by an axially symmetric cigar-shaped laser trap, elon- 
gated along the z-axis. The resulting trapping potential 
is U(r, z) = 0.5™[u; 2 r 2 + uj\z\ with uj r = 2ir x 437 Hz 
and uj z = 2ir x 27.7 Hz. The trapped Fermi cloud is ini- 
tially bisected by a blue-detuned beam which provides a 
repulsive knife-shaped potential. This potential is then 
suddenly removed, allowing for the two separated part 
of the cloud to collide with each other. The system is 
let to evolve for a given hold time t, then the trap is re- 
moved in the radial direction, and the system is allowed 
to evolve for another 1.5 ms during which the gas expands 
in the r-direction (during this extra expansion time, the 
confining trap frequency along the z-axis is changed to 
bj z = 2tt x 20.4 Hz), and finally a (destructive) image 
of the cloud is taken. The process is repeated from the 
beginning for another different value for the hold time 
t. The experimental results are eventually plotted as ID 
integrated density profiles for the different hold times in- 
vestigated: in Fig. 2 of Ref. [8|, the experimental ID 
density profiles at different times t are shown along the 
long trap axis. 

We simulated the whole procedure within the frame- 
work discussed above. As in Ref. §, we choose the 
initial density profile in the form of a static solution of 
hydrodynamic equations: 

n(r , M = 0> = fi(1 -^-|-^V (6) 

where h = [(2m^G/fi, 2 )/£] 3 / 2 /(37r 2 ). V rep represents 
an optically generated knife-shaped repulsive potential 
used to initially split the Fermi cloud into two spa- 
tially separated components that are led to collide with 
each other upon removal of such potential. V re p = 
V exp(— (z— z ) 2 /cr 2 ), where Vq — 12.7 [iK, a z = 21.1 fim 
and zq = 5/im. In Eq|5] we use the same values as in 
Ref.[8|, which provide a fit to the observed initial ex- 
perimental cloud density profile immediately after the 
removal of the knife potential. Here R z = 220 fim and 
R = 14 fim. In particular, the chosen value for the chem- 
ical potential fie — 0.53 fiK corresponds to a total of 
N = 2 x 10 5 6 Li atoms. 

During the time evolution of our system, when the two 
clouds start to overlap, many ripples whose wavelength 
is comparable to the interparticle distance are produced 
in the region of overlapping densities. These ripples also 
propagate backwards towards the trap boundaries, affect- 
ing larger and larger portion of the simulated cloud. We 
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FIG. 1: Density profile after t=1.5ms. The two curves show 
the smoothed and non-smoothed calculated profiles, respec- 
tively. The normalized density is in units of 10~ 2 /fim per 
particle. 



stress that these effects are quite similar to those found 
in a recent experiment Q by merging and splitting Bose- 
Einstein condensates. 

In order to properly compare our results with the ex- 
perimental data of resonant fermions Q , which are char- 
acterized by a finite spatial resolution, we smooth the 
calculated profiles at each time t by local averaging the 
density within a space window of ±5 /jm centered around 
the calculated point. This procedure will give smoothed 
density profiles with a spatial resolution close to the one 
characterizing the experimental setup of Ref . Q . 

The effect of smoothing is shown in Fig[TJ where the 
simulated density profile during the time evolution of the 
colliding clouds is shown before and after the local aver- 
aging procedure is applied. We wish to stress that this 
smoothing procedure is just a post-processing of the data 
obtained by the time evolution of the NLSE correspond- 
ing to Eqs. (QJ and (J5J) and therefore it does not affect 
the time evolution itself. 

The results of our simulations (solid lines), for the 
whole time duration of the experiments, and after the 
smoothing procedure is applied to the (y-averaged) den- 
sity profile at each time, are shown in Fig. [2] plotted 
along the long trap axis, for the same time frames as in 
the experiment. The experimental results (dotted lines) 
are also shown for comparison. 
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FIG. 2: ID density profiles at different times t showing the 
collision of two strongly interacting Fermi clouds. Solid lines: 
our calculations with no adjustable parameter. Dotted lines: 
experimental data from Ref. The normalized density is 
in units of 10 -2 /(im per particle. 



The time value shown in each frame corresponds to the 
time evolution of the initial profile, Eq. ([S]), before the 
trap in the radial direction is removed to let the system 
evolve for another 1.5 ms in the axial trap only, as done 
in the experiment. 

Note the striking correspondence between the experi- 
mental data and our simulation. We emphasize the fact 
that our simulations do not have adjustable parameter 
to be used to fit the experimental data, at variance with 
the model calculations presented in Ref. Q. Even after 
the smoothing procedures is applied, our simulated den- 
sity profiles exibit short-scale structures superimposed to 
the body of the cloud profiles. Such ripples are indeed 
present also in the experimental data, whereas in the 
model profiles used in Ref. [§] to fit the observed im- 
ages such oscillations are completely absent, due to the 
presence of a strong viscosity term in their model. 

The numerical results shown in Fig. [5] have been ob- 
tained by using A = 1/4, as previously discussed. In or- 
der to check the robustness of our results upon a different 
choice for A, we have performed various time-dependent 
calculations with the same initial conditions but using 
different values for A. It turned out that changing A from 
the optimal value A = 1/4 has profound consequencies on 
the long-time evolution of the colliding clouds, providing 
density profiles which are completely different from the 
experimental ones. 

This finding is not simply a numerical result, but has 
an important physical bearing. In fact, a strong depen- 
dence on A of the time evolution of a Fermi cloud made of 
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FIG. 3: "Soft" collision. The initial density profile is shown 
with dashed lines. 

a large number of atoms is at first sight surprising, given 
the fact that the gradient term should become less and 
less important with increasing N. We believe that such 
dependence is due to the presence of shock waves (i.e. 
regions characterized by large density gradients) in the 
colliding clouds. To check this hypothesis one should per- 
form experiments with different initial conditions, which 
do not evolve into shock waves. As an example, we have 
simulated a "soft" collision, where shock waves are not 
expected, by considering a system with a smaller den- 
sity of fermions (N = 4000 in the same trap used in the 
simulations of Fig. [5]) and where the two initial clouds 
are largely overlapping at the beginning of the simula- 
tion (this could be obtained experimentally by reducing 
the heigth/width of the "knife potential"), thus reducing 
the velocity of the impinging clouds. It turns out that 
the long-time behavior of such system is indeed indepen- 
dent on the chosen value of A. However, even for such a 
system, shock waves will eventually occur. As expected, 
such waves break into ripples whose wavelength depends 
on A. This is illustrated in Fig. [3l where the initial 
overlapping clouds are shown, together with the profile 
after t = 10 ms. For t < t s ~ 9 ms the shape of the 
time-evolved cloud is exactly the same for all values of A. 
After the shock time t s , however, the density profile de- 
velops a steep density gradient at the cloud boundaries, 
visible in the figure, which eventually breaks into a train 
of ripples. The three images shown in Fig[3] are taken 
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FIG. 4: Ripples wavelength I as a function of gradient co- 
efficient A of the density functional. The solid line shows the 
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just after the ripples have been produced. 

These ripples are characterized by a well-defined wave- 
length Z, which we estimated from the computed profiles. 
We plot in Fig. 0] the calculated ripples wavelength I 
as a function of A. We note that the calculated values 
for I follows very closely a y/\ law, which is expected on 
the basis of dimensional analysis by balancing the energy 
associated to dispersion to that of the nonlinearity. 

In conclusion, we have numerically studied the long- 
time dynamics of shock waves in the ultracold unitary 
Fermi gas. We have described the system by using an 
extended density functional approach, which has been 
used recently to successfully describe a number of static 
and dynamical properties of cold Fermi gases. Two main 
results emerge from our calculations: a) at zero tem- 
perature the simplest regularization process of the shock 
is purely dispersive, mediated by the quantum gradient 
term, which is one of the ingredient in our DF approach; 
b) the quantum gradient term plays an important role 
not only in determining the static density profile of small 
systems, where surface effects are important, but also in 
the fast dynamics of large systems, where large density 
gradients may arise. 

Finally, we stress that dispersive shock waves with a 
characteristic wavelength should be observable, accord- 
ing to our simulations, by using a soft-collision setup. 

We thank James Joseph and John E. Thomas for useful 
comments and suggestions, and for having sent us their 
experimental data. 
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